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Abstract 

We investigate the equal-mass 3-body system in general relativistic lineal gravity in 
the presence of a cosmological constant A. The cosmological vacuum energy introduces 
features that do not have a non-relativistic counterpart, inducing a competing expan- 
sion/contraction of spacetime that competes with the gravitational self-attraction of 
the bodies. We derive a canonical expression for the Hamiltonian of the system and 
discuss the numerical solution of the resulting equations of motion. As for the system 
with A = 0, we find that the structure of the phase space yields a rich variety of in- 
teresting dynamics that can be divided into three distinct regions: annulus, pretzel, 
and chaotic; the first two being regions of quasi-periodicity while the latter is a region 
of chaos. However unlike the A = case, we find that a negative cosmological con- 
stant considerably diminishes the amount of chaos in the system, even beyond that 
of the A = non-relativistic system. By contrast, a positive cosmological constant 
considerably enhances the amount of chaos, typically leading to KAM breakdown. 
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1 INTRODUCTION 



The N-body problem - that of describing the motion of a system of N particles interacting 
through specified forces - is one of the oldest problems in physics. Even now it continues to 
be a problem in fields ranging from nuclear physics to stellar evolution and cosmology. When 
the problem is that of N particles interacting through their mutual gravitational attraction 
it is particularly challenging. In the Newtonian case in three spatial dimensions, an exact 
solution is only known for the N=2 case, and for the general relativistic case in three spatial 
dimensions there is no exact solution for the equations of motion for even the N=2 case 
(although approximation techniques exist [1]). This is largly due to the dissipation of energy 
in the form of gravitational radiation, which results in the system needing to be solved by 
approximate solutions. 

Considerable progress has been made in recent years by studying systems with reduced 
spatial dimensions. Nonrelativistic one-dimensional self-gravitating systems (OGS) of N par- 
ticles have been very important in the fields of astrophysics and cosmology for over 30 years 
[2]. Although they are primarily used as prototypes for studying gravity in higher dimen- 
sions, there are physical systems with dynamics closely approximated by the one dimensional 
system. Very long-lived core-halo configurations, reminiscent of structures observed in glob- 
ular clusters, are known to exist in the OGS phase space [3]. These model a dense massive 
core in near-equilibrium, surrounded by a halo of high kinetic energy stars that interact only 
weakly with the core. Also, the collisions of flat parallel domain walls moving in a direction 
perpendicular to their surfaces and the dynamics of stars in a direction orthogonal to the 
plane of a highly flattened galaxy are approximated by the OGS system. 

In this paper we continue an ongoing investigation into relativistic one-dimensional self- 
gravitating systems (ROGSs). This is carried out in the context of a (l+l)-dimensional 
theory of gravity (lineal gravity) that models (3+1) general relativity by setting the Ricci 
scalar R equal to the trace of the stress energy of prescribed matter fields and sources. For 
the ROGS these sources are point particles minimally coupled to gravity. Hence, as in (3+1) 
dimensions, the evolution of space-time curvature is governed by the matter distribution, 
which in turn is governed by the dynamics of space-time [4]. Referred to as R = T theory, it 
is a particular member of a class of dilaton gravity theories on a line. What makes this theory 
of particular interest is that it has a consistent nonrelativistic (c — > oo) limit [4] (a problem 
for generic (l+l)-dimensional gravity [5]), and in this limit reduces to the Newtonian OGS 
as a special case. Furthermore, when the stress energy is that of a cosmological constant, it 
reduces to Jackiw-Teitelboim theory [6]. 

The three body ROGS has been studied previously and compared to the corresponding 
non-relativistic system [7] in the case of a zero cosmological constant for both equal and 
unequal masses [8, 9]. The degrees of freedom in both the non-relativistic and relativistic 
systems can be rewritten in terms of a single particle moving in a two dimensional potential 
well of hexagonal symmetry. In the non-relativistic case the potential grows linearly as a 
function of radial distance, but in the relativistic case the growth is non-linear, yielding a 
potential in the shape of a hexagonal carafe. For both systems, two broad categories of 
periodic and quasi-periodic motion were found, referred to as annulus and pretzel orbits, as 
well as a set of chaotic motions appearing in the phase space between these two types. The 
phase space between the two systems was found to be qualitatively the same despite the 
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high degree of nonlinearity in the relativistic system. 

In this paper we study the 3-body ROGS with non-zero cosmological constant A. This 
situation has no non-relativistic analogue. The effect of the cosmological constant is to induce 
a competing expansion/contraction of spacetime that competes with the gravitational self- 
attraction of the bodies. This system has been studied in some detail in the case of two 
interacting bodies, and it has been shown that relativistic effects are considerably enhanced 
by the presence of A [10]. For example a sufficiently large and positive A can overcome the 
attraction of the bodies, which then "lose casual" contact over a finite amount of proper 
time. The effect of A on the chaotic behaviour of the 3-body system has never been studied. 
It is the purpose of this paper to examine this problem, restricting our investigation to the 
equal- mass case. 

In Sec. II we review the formalism of the N-body problem in lineal gravity. We discuss 
the canonical decomposition of the action and find the Hamiltonian of the particles in terms 
of the dilaton field, which is determined by a set of constraint equations. We employ two 
different methods for solving the constraint equations, which allow us to solve the equations 
of motion for systems with both positive and negative values for the cosmological constant 
in Sec. III. In Sec. IV we define a coordinate system that allows us to describe the three 
particle system as a single particle on a plane, moving in a potential well with hexagonal 
symmetry. We describe the shape of this potential and its dependence on momentum and 
the cosmological constant. Here we also find upper and lower bounds on allowed values of 
the cosmological constant in a system with a given energy. Using these new coordinates, 
we describe our methods for numerically solving the system in Sec. V, and our methods of 
obtaining orbits, Poincare maps, and graphs of the oscillation patterns of the three particles. 
In Sec. VI, we numerically solve the equations of motion in the equal mass case with 
different values for the cosmological constant. We find two broad categories of periodic and 
quasiperiodic motions that we refer to as annulus and pretzel patterns, as well as a set of 
chaotic orbits found in the region of phase space in between these two types of orbits. We 
also identify general changes to the system in the presence of both a positive and a negative 
cosmological constant. Finally, in Sec. VII we present various Poincare maps and discuss 
how the global phase space is distorted in the presence of a non-zero cosmological constant. 
We find that a negative cosmological constant considerably reduces the amount of chaotic 
behaviour, whereas a positive cosmological constant considerably enhances it. In Sec. VIII 
we discuss the salient features of our solutions and make some conjectures regarding their 
general properties. We close the paper with some concluding remarks and directions for 
further work. 
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2 Canonical Reduction of the N-body Problem in Lin- 
eal Gravity 

As in previous work [7, 8, 9] we begin with the action integral for the gravitational field 
coupled with N point particles, which is 

r T 1 r 1 ~\ N r ( rlz^ dz v 1 5 

I = Jd 2 x _^|*i2+-^v#.*V^ + A}-E^«y ^«{-^)^r^r| S 2 (x-z a (r a )) 

(1) 

where g^ v and g are the metric and its determinant, R is the Ricci scalar, r a is the proper time 
of the a-th particle, k = 87rG/c 4 is the gravitational coupling, and with a scalar (dilaton) 
field This action describes a generally covariant self-gravitating system (without collision 
terms, so that the bodies pass through each other), in which the scalar curvature is sourced 
by the point particles and the cosmological constant A. Variation of the action with respect 
to the metric, particle coordinates, and dilaton field yields the field equations 

R — A — kT^ (2) 
d i dz v a 1 v ( \ d z a dz% 

d^W a } + T ^ {Za) ^~ a d^ a =0 (3) 

\ Vv^ ~ V A * Va * - V 2 ^) - V^V^ = + 2 9tw (4) 
where the stress-energy due to the point masses is 

N f 1 dz a dz p 

T lv = H m a \ dT a ^=g (ia g vp a a 5 2 (x - z a (r )) (5) 
0= i J V~9 dr a dT a 

and is conserved. We observe that (2) and (3) form a closed system of N+l equations for 
which one can solve for the single metric degree of freedom and the N degrees of freedom of 
the point masses. The evolution of the dilaton field is governed by the evolution of the point 
masses via (4). The left-hand side of (4) is divergenceless (consistent with the conservation 
of T^ u ), yielding only one independent equation to determine the single degree of freedom of 
the dilaton. 

We make use of the decomposition y^gR = —2d (^K) + 2di(^/ 7 yN 1 K — 7 _1 <9iiVo) 
where the extrinsic curvature K = (2V 7)~ 1 (2(9iV 1 — 7~ 1 iV 1 <9i7 — <9 7), and rewrite the 
action in the form 



1 = Jdx 2 {X>^ (x - z a (x )) + 7T 7 + IMf + N R° + A^i? 1 j 



(6) 



where 7 = gu, N = (—g 00 ) 2, N± — g w , ir and II are conjugate momenta to 7 and \& 
respectively. The quantities N and Aq are Lagrange multipliers that enforce the constraints 
BP = = R\ where 



fl » = - KV? V + 2 K ^n + gL-(-^j + ^-^ + ml s( x - Za { x ")) ( 7 ) 
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V^-^' + M + ^ix-z^)) (8) 

with the symbols ( ' ) and ( ' ) denoting d and di, respectively. This action leads to the 
following system of field equations: 



n + N { — V7vr 2 - — yrn + (tff - -A_ - V f a 5 (x - z a (t)) 



+ Nl \-- 2 U^' + - + Y,- 2 5(x-z a (*))} + K-^—V + N[— = (9) 

{ 7 2 7 a 7 2 J 2/t^77 7 



7 - iV (2/t v ^ 7 7T - +iVi— -2iV( = (10) 

7 



= (11) 



R 1 = (12) 



n + ^ ( --Aqn + — !-= W + ) = o (13) 

V 7 2 Kv /7 /t^/7 y 



* + ^o(2kv^7t) -JVi =0 (14) 



2 M + m 2 7 7 7 dZa 

V 7 a 



— iVi 

z a -N I +^ = (16) 

/^ +m 2 7 

We can solve the constraint equations (11) and (12) in terms of the quantities / \f~j)' 
and 7r', since they are the only linear terms present. The generator obtained from the end 
point variation can then be transformed to fix the coordinate conditions. We can consistently 
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choose the coordinate conditions 7 = 1 and II = [11]. Eliminating the constraints, the 
action then reduces to 

1 = Jd 2 x ^p a z a 5 (x - z a ) - W j (17) 



where the reduced Hamiltonian is 

H 



J dxH = -i J dxA^ (18) 



where A = d 2 /dx 2 , and * = * ( x, z a ,p a ) and is understood to be determined from the 
constraint equations which are now 

A*-^ + M-± + K^ti + ml6(x-z a ) = (19) 
2A X + y £ l PaS(x-z a ) = (20) 

a 

where ir = X ' ■ The consistency of this canonical reduction can be demonstrated by show- 
ing that the canonical equations of motion derived from the reduced Hamiltonian (18) are 
identical with (15) and (16) [11]. 



3 Solving the Constraint Equations 

We can solve these two equations exactly for iV = 3 particles by solving the equations in 
the empty spaces between particles and then matching the fields at each of the particle 
positions, imposing the condition that the Hamiltonian remains finite as 1 — > 00. After 
these constraints are applied we are left with an equation containing only one unknown. 

We divide the space into four regions determined by the arbitrary particle positions 
chosen so that z 1 < z 2 < z 3 . We label the four regions: x < Z\ ((1) region), Z\ < x < z 2 , 
((2) region), z 2 < x < z 3 ((3) region), and z 3 < x ((4) region). In each of these regions the 
solution to (20) is simply 

1 3 

X= --^2p a \x-z a \-eXx + eC x (21) 

4 a=l 

with X and C x being two integration constants. The factor e (e 2 = 1) changes sign under 
time reversal; it has been introduced to explicitly manifest the property that x a l so changes 
sign under time reversal. We find an easier expression for (19) by using the substitution 

^ = -4 log |0| (22) 

which gives us 

A0 - \<t> {~ + (X'f] = 1 E Vrf + <4> (*a) S(X - Za) (23) 
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In the four regions defined above equation (23) has the solution 



Aie? KlX + B\e i Klx in the (1) region 

A 2 e^ K2X + B 2 e~^ K2X in the (2) region 

A 3 e^ KzX + B 3 e~i K3X in the (3) region 

A^ Kax + B±e~i KAX in the (4) region 



(24) 



where the K; are defined as 



K 2 = 



jn 2 


X ~ l( P l + P 2 + P3) 


2 A 
2 


K 2 ' 


X - i(~Pi + p 2 + Pa) 


I 2 A 
2 



X + f(Pl +P2 -Pz) 



A 
2 



X + - 4 (P1+P2+P3) 



(25) 



The matching conditions that must be satisfied require that is continuous and its derivatives 
are consistent with integration of eq. (23) across the location of each particle: 

(26) 





02 (4 ) 


= 0i (*r) 




03 (4) 


= 02 


02 (4) 


04 (4) 

- 0i (*r) 


= 03 (*a) 

= ^\/pi + m i0( z i) 


03 (4) 

01 (4) 


- 02 (W) 

- 03 ( 2 3~) 


= ^VP2 + ^20(^2) 
= ^ \/P3 + ™30 (^) 



where a (zj) = lirn^^^ a (x). 

This leaves us with six equations for ten unknowns (A% S's, X, and C x ). There is also 
the condition that the Hamiltonian remains finite, which can be done by requiring 

^ 2 - 4«V + 2Ax 2 = C±x (27) 

where C± are constants to be determined [11, 12]. This gives 

A 1 = B A = 

and four other equations that contain C±. This leaves us with ten equations for ten unknowns 
(A 2) A 3: At, Bi, B 2 , B 3 , C x , C_, C+, X). At this point, one can do the integration explicitly. 
The Hamiltonian becomes 



H = — 



K 



dxA^J = -- 



1 



K 



4(f)' 

'T 



+oo 



2 (K, + K A ) 



K 
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Without loss of generality, we choose the centre of inertia frame Pi + p 2 + Ps = 0, which 
simplifies the Hamiltonian to 



A 



K 



(28) 



The Hamiltonian can only depend upon the relative separation of the particles and their 
conjugate momenta. Furthermore, in the equal mass case, the system is symmetric under 
particle interchange. We introduce the following notation 



Zij — (Zi Zj) 

Sij = sgn (zij) 



R 



i± 



\ 



X ~ 7 [12 P aSai ± Pi 
4 \a=l 



A 
2 



L* = 



M tj = K^Jp 2 , + mf + 2 Sij (R i+ - Ri_) 
Li = -n^p 2 + m 2 t + 2 (R i+ + /?,_) 

s iJ s ik ^pf + m 2 + 2 (R i+ + Ri_) 

j<k^i J 



(29) 
(30) 

(31) 

(32) 
(33) 

(34) 



where these quantities are defined in such a way that they automatically take care of the 
crossings of particles via s^. Using them, the derivation of the full determining equation for 
the Hamiltonian is similar to that for the A = case [8] ; the result is 



LiL 2 £ 3 = M 12 M 2 iL*e^ 2l[(Ml2+Ll)231 - {M21+L2)z321 
+M23M3 2 Lie^ 32[(M23+L2)2l2_(M32+L3)zi3] 
+M3iMi 3 L2e^ Sl3 ^ M31+L3 ^ 23 ~ ( ' Ml3+Ll ' )Z2 ^ 



(35) 



or more compactly 



LiL 2 L 3 = J2 leij-fclA/ij-MjiL^e^yl^+^^i-Wi+^^w] 

ijk 



(36) 



where is the Levi-Civita tensor. 

While technically correct, the determining equation (35) is not useful computationally 
insofar as some of its terms can become imaginary with a sufficiently large cosmological 
constant; this is a consequence of eq. (25). However we can rewrite the determining equation 
to circumvent this problem. Consider the case where z\ < z 2 < z 3 . Redefining the following 
terms 



M i = K-Jpf + m 2 



(37) 



K 3 = 



1 



-, 2 



X + 4 [12 a 3iPi 



A 
2" 
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where = — 1 when j < i and = 1 when % < j we can get a new form for the determining 
equation 

((Mi + ^i) (m 3 + K 4 ) M 2 + (Mi + £1) iT 3 2 + (M 3 + £ 4 ) K 2 2 ) tanh (^3^32) tanh (jK 2 z 2 i 



+ ((Mi + M 2 + Xi) (M 3 + £4) + £T 3 2 ) K 2 tanh Qtf 3 *32 
+ ((Mi + k t ) (M 2 + M 3 + £ 4 ) + K 2 ) X 3 tanh Qtf 2 z 21 

+ (Mi + m 2 + m 3 + + k 4 ) k 2 k 3 



= (38) 

When the center of momentum is set to zero only the terms K 2 or K3 can be imaginary. 
Consequently (38) is either purely real, or purely imaginary (in which case the % can be 
factored out). Permutation of the particles results in the exact same determining equation 
with the indicies appropriately switched. 

The Hamiltonian H is only implicitly determined from (38). Fortunately we do not need 
its explicit form since we can implicitly take derivatives on both sides of the determining 
equation (38) with respect to the phase space variables and from there extract the canonical 
equations of motion. From (28) with p z = ELi Pa = 0, we obtain 



^ = t A. K 2 X 2 _ A = dX 
dx k K\dx k y 2) ^tfx 2 - f 9x k 

where x k is any of the canonical variables. We then have the simple but tedious task of 
taking ^ of both sides of equation (38), collecting the derivatives of X, and then converting 
them to derivatives of H via (39). This yields the canonical equations of motion 

dH , . 

Za = 7T- 40 
dp a 

Pa = ~ (41) 

dz a 

where the dot denotes a derivative with respect to the coordinate time. These equations of 
motion are straightforward to calculate but highly tedious and will not be included here. 
When calculating the equations of motion for the particles when they are not in the arrange- 
ment Z\ < z 2 < z 3 we temporarily change the labels of the particles so that they do satisfy 
this condition. This allows us to use the derivatives of (38) to calculate the equations of 
motion for each of the particles, after which the original particle labels are returned. This 
method works regardless of whether or not the particles are identical. 



4 General Properties of the Equations of Motion 

We can now numerically calculate the equations of motion for each of the three particles. 
This would give us six equations of motion when our physical system only has four effective 
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degrees of freedom: two relative separations and their two conjugate momenta. Introducing 
the new variables 



Z2 ~ Zi 



V2p (42) 



{zi + z 2 ) + 2^3 = v^A 
z\ + z 2 + z 3 = Z 



and their conjugate momenta 



Pp = ^(P2-Pi) (43) 
Px = (- (Pi + P2) + 2p 3 ) 
Pz = 7^ (Pi + P2 + Ps) 

which have been chosen to satisfy the canonical commutation relations {qi,Pj} = Sij. This 
choice of variables gives the Hamiltonian an explicit sixfold symmetry. We can set the 
conjugate variables Z and pz to arbitrary values since Z is irrelevant in our equations and 
we can fix the center of inertia by setting p z = without loss of generality. Inverting these 
relations, we get 

z 2 -z 1 = V2p (44) 

Z3-Z2 = -^(Vsx-p) 

Z3-Z1 = -^=(V3\ + p) 

and 




PI = "T^-Tr (45) 

P3 = 

This choice of variables allows us to write the determining equation in terms of only 
{p, \,p p ,px} for a given cosmological constant. In this way, the relativistic Hamiltonian can 
be regarded as a function H = H (p, \,p p ,p\) 

We see from this perspective that the linear three-particle system is equivalent to a single 
particle moving in a two dimensional "potential well". This allows the variables p,\,p p ,p\ 
to be interpreted as the coordinates and the conjugate momenta of this single particle, which 
we call the hex-particle due to the hexagonal symmetry of the potential. In the Newtonian 
case, the potential takes on the shape of a hexagonal-shaped cone with planar sides and is 
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independent of the momenta [13]. The relativistic potential well is obtained by regarding the 
potential to be the difference between the Hamiltonian and the relativistic kinetic energy 

V (p, A) \ Pp=a ,p x =b = H (p, A, p p = a,p x = b)- \J (mc 2 ) 2 + (c \p\) 2 (46) 

and is dependent on the momentum \p\ = \/a 2 + b 2 of the hex-particle as well as its position 
in the (p, A) plane 4 . 

We get our first look at the relativistic potential by considering the case when p p = = p\. 

V (p, A) \ Pp =o, Px =o = H (p, A, p p = 0, p A = 0) - mc 2 

At very low energies this relativistic potential is indistinguishable from the potential for the 
Newtonian case. However, as shown in table (1), even at energies only moderately larger 
than the rest mass, the sides of the hexagon become convex in the relativistic potential. 
Though it retains its hexagonal symmetry, the growth of the relativistic potential is very 
rapid as p and A increase. The size of the cross section of the potential reaches a maximum 
at an energy V Rc just over twice the rest energy of the particle, after which the diameter of 
the potential decreases like In (V R ) /V R with increasing V R [8]. 

The part of the potential on the branch with V R > V Rc is in an intrinsically nonperturba- 
tive relativistic regime. The motion for values of V R larger than this cannot be understood as 
a perturbation from some classical limit of the motion. The nonrelativistic hexagonal cone 
becomes a hexagonal carafe in the relativistic case, with a neck that narrows as Vr increases. 

When p p = = p\ it is straightforward to show that the potential is independent of the 
cosmological constant, a feature noted previously in the 2-body system [10, 12]. However, 
for nonzero p this is not the case, and we map out the potential to show how it changes as 
A varies. It is convenient to instead map the potential with set values of radial momentum 
p r and angular momentum pg of the hex-particle where 

Pr = P P cos 9 + pa sin 9 



Pe = —p P sin 9 + p x cos 9 



9 = arctan 

since the change to polar coordinates manifestly retains the hexagonal symmetry of the 
potential, which greatly simplifies the subsequent analysis. 

First, when A = and the hex particle has positive radial momentum p r , the width of 
the potential at lower energies increases with increasing p r and the sides of the hexagonal 

A Note that earlier definitions of the potential [8] defined it to be the value of the Hamiltonian at zero 
momenta. 
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A = Q Pr = Q,3 



A = Q Pr = Q 



A = p r = -0.3 

Table 1: Potential plot for A = from both side and top views with solid lines denoting 
equipotentials. 
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A = Q Pe = Q,5 



Table 2: Potential plot with positive angular momentum from both side and top views with 
solid lines denoting equipotentials. 



A = Q,l Pr = Q,3 



A = 0.1 p r = -0.3 

Table 3: Potential plots for A > from both side and top views with solid lines denoting 
equipotentials. 
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A = -0.1 p r = 0.3 



A =-0.1 p r = -0.3 

Table 4: Potential plots for A < from both side and top views with solid lines denoting 
equipotentials. 
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cross section become more convex at lower energies, resulting in a star shaped cross section, 
as shown in table (1). This widens the bottom of the potential, and lowers the value of 
Vrc at which the cross-section of V is at its largest. When the hex-particle has negative p r 
the opposite happens: the width of the potential decreases and the sides of the hexagonal 
cross section become more concave, resulting in a more rounded cross section. Though its 
sides are initially much steeper, the potential curves back in on itself much more gradually 
resulting in a higher value of Vr c . 

A similar analysis for nonzero angular momentum p e , as shown in table (2), indicates 
that only for very high values does it have any substantive effect on the potential . As pg 
increases the sides of the hexagonal cross section become concave and slightly skewed in the 
direction opposing the angular momentum. We expect that this is the cause of the previously 
observed rotation of the annulus orbits at higher energies [8] for the A = case. 

When A > the potential, shown in table (3), shifts down slightly when p r is positive. 
At lower energies the sides are pushed farther out while becoming more convex in a manner 
similar to that which occurs for positive radial momentum in the original (A = 0) potential. 
A side effect of this is that it seems to lower Vr c , which may cause the system to become 
unstable at lower energies. Associated with this is the appearance of a positive critical value 
of the cosmo logical constant, which we will discuss in more detail later. For negative p r and 
A > 0, at lower energies the potential is pulled inward and becomes more concave, again 
in a manner similar to the negative radial momentum in the A = case. However, the 
magnitude of this effect of A on p r < is much less than the effect on p r > 0. The effects 
noted above for angular momentum also seem to increase when A > 0. Overall, the presence 
of a positive cosmological constant seems to enhance the effects that momentum has on the 
A = potential. 

When A < 0, as shown in table (4), the opposite effect happens. For p r > the sides are 
pulled in relative to the effects of momentum in the A = potential, reducing its diameter as 
well as making the sides less convex. For p r < 0, the sides of the potential are very slightly 
pushed back out relative to the effects of momentum in the A = potential, increasing its 
diameter. The effects of angular momentum are also counteracted when A < 0, so overall, 
the presence of a negative cosmological constant seems to reduce momentum-dependence in 
the potential. 

We have found both upper and lower limits on the range of allowed values of the cosmo- 
logical constant A that are dependent on the energy in the system. As described before, in 
order to solve the equations of motion we must calculate the derivatives of X from Eq. (21) 
with respect to each Zi and pi, and then use Eqs. (28) and (39) to convert X-derivatives to 
iJ-derivatives, yielding the canonical equations of motion. However, both X and H must be 
real. Thus from eq. (28) we find that for any given value of H, A must be chosen in a range 
that satisfies 

H > -J 

"4 2 

resulting in a negative critical value for A 

H 2 k 2 

■''-negcrit g 

As A — > An egcr i t , we see that X — > . From eq. (39) all of the derivatives of H tend to 
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(47) 



zero, rendering the particle motionless. If we take a fixed value for H and increase it by an 
arbitrarily small amount h, when A is the negative critical value for an energy of H Q , from 
eqs. (28, 47) we get 



When h is very small, so is X, and we can take the approximation h ~ 8^- , yielding 



dH dh 16X8X 32hdX 



dx n dx n H dx n V dx r , 



(48) 



for the equations of motion, where x n is any of the canonical variables. For small h the 
equations of motion are all scaled by \fh, which is equivalent to scaling the time factor 
by -j=. Consequently, when the cosmo logical constant is very close to its negative critical 
value for a system with a given energy, the hex-particle will not so much be restricted in 
its movement in the (p, A) plane, and follow that orbit at a much slower pace. However, as 
we will describe in more detail later, we have found that for small negative values for the 
cosmological constant the frequency of the hex-particle's movement increases, so this trend 
only appears for very small h. 

We have not been able to identify analytically the positive critical value of A . Rather we 
identify it numerically: when running simulations with A too large the three particle system 
breaks out of its bound state, and our simulation breaks down. We have not yet been able to 
conclusively find why the system fails where it does, but we believe it may be due to the fact 
that, when radial momentum is present in the hex-particle system, a positive cosmological 
constant seems to lower the value of Vr c , rendering the system intrinsically nonperturbative 
where it was previously in a perturbative regime. This hypothesis is supported by the fact 
that, for the energies that we have tested, the maximum value of A allowed seems to decrease 
substantially as H increases. 

The effect on the hex-particle's orbits from approaching these positive and negative crit- 
ical values of the cosmological constant will be further discussed in later sections. 



5 Methods for Solving the Equations of Motion 

We find it easier to study the equal mass three body system by analyzing the motion of 
the hex-particle in the (p, A) plane. In this system the bisectors joining opposite vertices 
of the potential well correspond to two of the three particles crossing one another, causing 
a discontinuous change in the hex-particle's acceleration. Therefore, in the Hamiltonian 
formalism, the motion of the hex-particle is governed by a pair of differential equations which 
are continuous everywhere except at the three hexagonal bisectors p = 0, p — v^A = 0, and 
p + y^A = 0. These correspond to the crossings of particles 1 and 2, 2 and 3, or 1 and 3 
respectively. 

We allow the hex-particle to freely cross over each of these bisectors, which corresponds 
to a pair of particles passing through each other. An analogous system in the Newtonian 
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case, consists of the motions of a ball under a constant gravitational force elastically colliding 
with a wedge [14]. This system is nearly identical to the one we study insofar as, in the equal 
mass case, an elastic collision between a pair of particles in the three body system cannot be 
distinguished from a crossing of two equal mass particles. However, a distinction between 
the two systems in a certain class of orbits has been previously observed [8]. 

The nonsmoothness of the potential gives rise to interesting dynamics in the system. We 
shall refer to two distinct types of motion in our subsequent analysis [14]. A motion, corre- 
sponding to the same pair of particles crossing twice in a row in the three body system (or 
equivalently the hex-particle crossing a single bisector twice in succession), and B motion, 
in which a single particle crosses both of the other particles in succession (or equivalently 
the hex-particle crossing two successive hexagonal bisectors). For a given orbit we can char- 
acterize the motion by a sequence of these symbols with a finite exponent n denoting n 
repeats and an overbar denoting an infinite repeated sequence. It is important to note that 
the classification of a crossing as A or B motion is dependent on the previous crossing, so 
there is an ambiguity in the classification of either the final or the initial crossing. This is 
resolved by taking the initial crossing as unlabeled; since we are considering large sequences 
of motion this ambiguity causes no practical difficulties [8]. 

There is a technical difficulty in comparing trajectories in systems with different values 
for the cosmological constant. Such changes induce a change in the determining equation and 
so we cannot use exactly the same initial conditions. We will deal with this by comparing 
trajectories with the same fixed-energy (FE) initial conditions. This is done by fixing the 
initial values of H, p, A and p p , and adjusting p x so that the Hamiltonian constraint is 
satisfied. 

We have implemented three methods of analysis to study the motion of the system. 
First, we plot the trajectory of the hex-particle in the (p, A) plane and compare the motion 
under different initial conditions and in systems with a different value for the cosmological 
constant. Second, we plot the motion of the three particles as a function of time which gives 
us a different method for visualization, allowing us to identify differences in the various 
types of motion. Third, we construct Poincare sections by plotting the radial momentum 
(p r , labeled gainst the square of the angular momentum (pg, labeled as z) of the 

hex-particle each time it crosses one of the bisectors. When all particles have equal mass, all 
of the bisectors are equivalent so all of the crossings can be plotted on the same surface of 
section. This allows us to easily identify regions of periodicity, quasiperiodicity, and chaos, 
and we will discuss each of these features. 

There is no closed-form solution for either the determining equation (35) or to the equa- 
tions of motion so we must solve these equations numerically. We proceed by introducing 
dimensionless variables 

X = M tot x H = M tot h A = M tot KA 
where M tot = J2a=i m a, and also 

4 

Pi = M tot pi Zi = -^ri—Zi 
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Substituting these variables into Eqs. (31), (32), (33), (34) we get 



Ri± = kM ( 



tot ^ • 




\fpl + m l = M tot\jv\ + m 2 a 



^ = KM tot R i± 



which gives us 



Mi,- = KM tnf M iA U = KM tot Li L* = KM tnt L 



Hot J 



The notation for these new variables have exactly the same form but without the factor of 
kappa inside the square roots and all multiplied by a factor of kM w . Thus when writing the 
determining equation in dimensionless variables the factors of («M tot ) 3 cancel out everywhere 
and the exponential factors become 



3 % 



KM tot (Mij + Li) -^j-z k i - KMua (Mji + LA 

V / KMua V J} KM tot 



(M i3 + z ki - (Mji + Lj 



So we see that introducing dimensionless variables is equivalent to writing the determining 
equation with K,M tot set equal to unity. 

We numerically integrate the equations of motion generated from the determining equa- 
tion (38) using a MatLab ODE routine (ODE45 or ODE113). When the particles are not in 
the arrangement z\ < z-i < z% our program temporarily changes the labels on the particles 
so they satisfy this condition. The equations of motion are then calculated for each of the 
particles and their original labels are returned. This method works regardless of whether or 
not the particles are identical. 

For the Poincare sections we stopped the integration each time the hex-particle crossed 
one of the bisectors by using an "events" function and saving the values of radial and angular 
momentum for plotting. Ideally for each chaotic trajectory the Poincare section should be 
allowed to run for a long time in order to accurately determine which regions of the plane 
the chaotic motion extended to and which regions were restricted. 

In order to control errors, we imposed absolute and relative error tolerances of 10~ 8 for 
the numerical ODE solvers. For the values of H we studied this yielded numerically stable 
solutions. We tested this by checking that the energy in the system remained constant 
throughout the motion to a value no larger than 1CT 6 . 



6 Equal Mass Trajectories 

We look at the cases where all three particles have equal mass and the hexagonal symmetry 
of the potential well is maintained. As described in previous papers [8] we find three general 
classes of orbits which we denote as annuli, pretzel, and chaotic. Examples of these classes 
of orbits have been found at all energies and in the presence of a positive and negative 
cosmological constant. 

Furthermore, in each class there are orbits that eventually densely cover the region of 
the (p, A) space that they occupy and orbits that do not. The latter situation corresponds 



17 



Figure 1: Examples of near-periodic orbits (run for 200 time units), trajectories do not cover 
the entire (p, A) space. Periodic annulus and pretzel orbits can be found for almost all 
testable energies and values of A. 

to a regular orbit that repeats itself after a finite time - the symbol sequence consists of a 
finite sequence repeated infinitely many times, resulting in a periodic trajectory. The former 
situation corresponds to an orbit where the same finite sequence is repeated, but with A 
motion added or removed at irregular intervals, resulting in a quasiperiodic trajectory. 

The quasiperiodic orbits closely resemble the periodic orbits, except that they fail to 
exactly repeat themselves and have some form of precession that causes the orbits to fill a 
region of space. Thus these quasiperiodic orbits show a high degree of regularity, manifest 
by its periodic symbol sequence. 

6.1 Annulus Orbits 

The annuli are orbits where the hex-particle never crosses the same bisector twice in succes- 
sion, resulting in the symbol sequence B and an orbit in the shape of an annulus encircling 
the origin in the (p, A) plane. 

These annulus orbits can be both quasiperiodic and fill in the entire ring, while a few 
repeat themselves after some number of rotations about the origin. Since periodic orbits are 
very difficult to find numerically, the orbits in the figures are actually just very close to being 
periodic, so that the pattern of the periodic orbit can be seen. 

There is a slight rotation in the annulus orbits as H increases, but otherwise little quali- 
tative difference between the annulus orbits occurs for different energies. However we do see 
some qualitative changes to the orbits as A changes. For the larger positive values of A that 



18 



Figure 2: Arbitrary examples of quasi-periodic orbits that densely fill the (p, A) space (run 
for 200 time steps). Quasi-periodic annulus and pretzel orbits can be found at all testable 
energies and values of A. 
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Figure 3: Annulus orbits for zero cosmological constant (H = 1.2 lower left; H = 1.8 upper 
right) shown in conjunction with their corresponding three-particle trajectories. The vertical 
axis is the displacement from the origin in units of K,M tot c 2 /A. These orbits have been run 
for 200 time steps but the three-particle trajectory plots have been truncated after 30 time 
steps. The rotation in the annulus orbits is greater for the higher energy example. 
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Figure 4: Overlay of annulus orbits at H = 1.5 starting with the same FE initial conditions 
but with different values for A run for 200 time steps and their corresponding three-particle 
trajectories truncated after 30 time steps (top A = 0.02; middle A = 0; bottom A = —0.25). 

we can obtain numerically, we find for the same FE initial conditions that the outer part 
of the annulus extends farther out into the (p, A) plane. We also find that the width of the 
annulus also increases more than if the orbit were simply linearly dilated (ie 'photographi- 
cally enlarged'). Conversely for large negative A the orbit covers less of the (p, A) plane and 
the annulus becomes thinner by a greater factor. These effects can be seen in figure (4) and 
seem to be consistent with how the potential changes with A. For positive A, the effect that 
positive radial momentum has on pushing out the potential and negative radial momentum 
on pulling the potential in towards the origin is exaggerated. Hence when the hex-particle is 
moving radially outward it will go farther because the potential is being pushed back. Simi- 
larly when the hex-particle is moving radially inward it will go farther because the potential 
is being pulled in and driving it farther. For negative A the opposite occurs and the radial 
motion is somewhat damped, resulting in a smaller range of radial motion. 

Another characteristic change induced by A is that the frequency of the three particle 
motion decreases as A becomes positive and increases as A becomes negative. At values 
of A extremely close to the negative critical value given in Eq. (47) the frequency of the 
motion should eventually decrease as described in Eq. (48). However, our simulation breaks 
down due to error tolerance limits so at this time we have not been able to verify this 
property numerically. This non-linear change in frequency can result in the same set of 
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Figure 5: Examples of pretzel orbits for A = (if = 1.8, lower left; H = 1.5, upper right) 
run for 200 time steps and shown in conjunction with their corresponding three-particle 
trajectories truncated after 60 time steps. Interactions resembling a quasistable bound 
subsystem of two particles orbiting a third particle can clearly be seen. 

initial conditions producing periodic or quasiperiodic orbits depending on A. This change in 
frequency has been observed in all three types of orbits. 

6.2 Pretzel Orbits 

In Pretzel orbits, the hex-particle essentially oscillates back and forth about one of the bisec- 
tors, corresponding to a stable or quasistable bound subsystem of two particles. This bound 
pair then orbits the remaining particle, which corresponds to the hex-particle oscillating 
along the same bisector. The existence of a two-particle bound subsystem was discovered 
in the Newtonian case [15] and has been observed previously in the A = case [8]. Sym- 
bolically, these orbits can be written as Yiijk (A ni B 3m :>) k , where ni,rrij,lk £ Z + , with some 
Ik possibly infinite. This motion results in an extremely diverse collection of orbits. Many 
families of regular orbits exist, containing one base element (for example AB 3 ) and a se- 
quence of elements formed by appending an A to each existing sequence of A's (for example 
{AB 3 , A 2 B 3 , A 3 B 3 , ...}). This results in an extremely complex structure for the phase space, 
which we will analyze later. 
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Figure 6: Overlay of pretzel orbits starting with the same FE initial conditions but with 
different values for A run for 200 time steps and their corresponding three-particle trajectories 
truncated after 30 time steps (top A = 0.02; middle A = 0; bottom A = —0.25). 

The B 3 term in the above sequences corresponds to a 180° swing of the hex-particle 
around the origin, and the orbits in the (p, A) plane that result from this motion comprise 
a wide variety of twisted, pretzel-like figures. Again, in this class of orbits we find both 
periodic orbits with an infinitely repeating symbol sequence, and quasiperiodic orbits that 
densely fill a cylindrical tube in the (p, A) plane which usually has a kink about the A = 
line. 

Qualitatively, the characteristics of the pretzel orbits are very similar at different energies. 
At all of the values of H that we can obtain numerically we can find orbits of all of the 
different allowed families that have regular symbol sequences. One of the differences is that 
it seems that the kinks in at the A = line become more pronounced at higher energies, as 
previously observed in [8]. 

As we change the cosmological constant we see many similar changes to the pretzel orbits 
that we found for the annulus orbits. As A increases the orbits extend farther out on the 
(p, A) plane, corresponding to the widening of the potential with positive radial momentum 
described before. Alternatively, as A decreases the orbits cover less of the (p, A) plane, 
corresponding to the contracting of the potential with negative radial momentum. We can 
find orbits with the same regularly repeating symbol sequence for all values of A that we can 
obtain numerically. However, as A becomes positive we have found that some pretzel orbits 
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begin to show a form of chaotic behaviour not seen in the A < cases. We will describe 
these cases in more detail later. 

Again in the three body figures we find that the frequency of the motion is increased and 
reduced for negative and positive A respectively. This non-linear change in frequency can 
cause the same set of initial conditions to result in either periodic or quasiperiodic motion 
depending on A. However in addition to this we find that the orbit shape and symbol sequence 
can be drastically changed as well. We will be able to see the overall structure of the change 
when we begin to look at the bottom region of the Poincare maps. 

6.3 Chaotic Orbits 

The chaotic orbits are those in which the hex-particle wanders between A motion and B 
motion in a seemingly irregular fashion. On the Poincare section these orbits appear as 
densely filled regions. These orbits eventually wander into all areas of the (p, A) plane 
allowed by the energy constraint, a trait not seen in either the annuli or pretzel orbits. 
Areas of chaos can be found in this system at the transition point between annulus and 
pretzel orbits, though the size of the chaotic region seems to be highly dependent on the 
cosmological constant, as will be shown later. 

Chaotic orbits have been found at all values of H that are numerically obtainable and 
have essentially the same basic properties. However these orbits are generally hard to find 
at any energy as they are very sensitive to initial conditions. The top left orbit of figure 7 
gives us an example of a chaotic orbit found in between the annulus and pretzel regions on 
the Poincare section. We can see that the orbit sometimes shows the characteristic motion 
of an annulus orbit and sometimes shows the motion of a pretzel orbit. The three-particle 
trajectories show us how the particles switch roles during their quasi-regular motion at 
irregular intervals resulting in the chaotic motion. 

As with the annulus and pretzel orbits, the chaotic orbits extend farther out into the 
(p, A) plane as A increases, and the frequency of the three particle motion increases. However, 
chaotic orbits are non-periodic by definition, so a change in frequency can never make this 
type of orbit periodic. 

The most significant change to the chaotic orbits at different values of A is that as A 
becomes negative, the initial conditions that resulted in a chaotic orbit will generally result 
in either an annulus or pretzel orbit. Conversely, we have found more initial conditions 
resulting in chaotic orbits when A is positive. Overall there seems to be a strong positive 
correlation between the value of the cosmological constant and the size of the subset of initial 
conditions that result in chaotic orbits. This trend is illustrated by figure (8) where for the 
same FE initial conditions are run with different values of A. The initial condition that 
results in chaotic orbits for A > instead result in a regular annulus orbit (in this case) for 
A < 0. 

In addition, we have also found separate regions of chaos within the pretzel region of the 
Poincare map when A is positive. For certain initial conditions resulting in pretzel orbits 
when A = 0, when we increase A we find that the hex-particle begins to irregularly transverse 
increasingly larger regions of the (p, A) plane. However, these chaotic pretzel orbits do not 
cover the entire (p, A) plane as do the chaotic orbits in between the annulus and pretzel 
regions. An example of this type of orbit is shown on the bottom left of figure (7). In these 
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Figure 7: Examples of chaotic orbits run for 200 time steps between the pretzel and annulus 
regions (top right) and within the pretzel region (bottom left) along with their corresponding 
three-particle trajectories truncated after 60 time steps. The sizes of the regions of chaos are 
highly dependent on the cosmological constant and the chaotic oribts in the pretzel regions 
are only found when A > 0. 
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Figure 8: Chaotic orbits for H — 1.2 starting with the same FE initial conditions but with 
different values for A run for 200 time steps (top left A = 0.06; top right A = 0; bottom 
A = —0.175). Not overlapped to emphasize change in amount of chaotic motion 
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chaotic pretzel orbits - unlike the previously described chaotic orbits where all three particles 
switch positions at irregular intervals - we can see from the three-particle motion that while 
two of the particles orbit each other closely in a bound subsystem, the third particle shows 
fairly regular motion. The chaotic motion of the hex-particle results from irregular motions 
in the two-particle subsystem. 

7 Poincare Plots 

We now consider the Poincare sections for this system. These are constructed by plotting 
the square of the angular momentum (p 2 6) labeled as z) of the hex-particle against its radial 
momentum (p r , labeled as x) each time it crosses one of the bisectors. Since we are looking 
only at the case when all three particles have the same mass, all of the bisectors are equivalent, 
and so we can plot all of the crossings on the same surface of section. This allows us to find 
regions of periodicity, quasiperiodicity, and chaos. 

Since this system is governed by a time-independent Hamiltonian with four degrees 
of freedom, the total energy is a constant of the motion and the phase space is a three- 
dimensional hypersurface in four dimensions. If an additional constant of motion exists the 
system is said to be integrable, and its trajectories are restricted to two-dimensional sur- 
faces in the available phase space. Since the trajectories can never intersect, that constraint 
imposes severe limitations on the types of motion that integrable systems can exhibit. Trajec- 
tories may be periodic (repeating themselves after a finite interval of time) or quasiper iodic. 
Since the trajectories are, by definition, comprised by the intersection of two two-dimensional 
surfaces, they will always appear as lines or dots for quasiperiodic and periodic orbits re- 
spectively, on the Poincare section. This is a sharp contrast to the case when all orbits can 
move freely in three dimensions in a completely nonintegrable system. The extra degree of 
freedom permits an orbit to visit all regions of phase space, and the system typically displays 
strongly chaotic behaviour. Such trajectories appear as filled-in areas on the Poincare map. 

When an integrable system is given a sufficiently small perturbation, most of its orbits 
remain confined to two-dimensional surfaces. However, it is possible that small areas of chaos 
will appear sandwiched between the remaining two-dimensional surfaces, which can grow as 
the magnitude of the perturbation is increased, eventually becoming connected areas on the 
Poincare section. This phenomenon is called a Kolmogorov-Arnold-Moser (KAM) transition 
[16]. Within these regions, islands of regularity may remain for quite some time and can 
have an intricate fractal structure [17], but the system will typically become almost fully 
ergodic for sufficiently large perturbations [18]. 

The outer bound of the Poincare section is limited by the energy of the system. The 
first noticable feature is that it is not symmetric about the p r = axis, but is instead 
skewed to the right, with the deformation increasing with H. This is at first puzzling; the 
trajectories of a subset of the annulus orbits always have positive radial velocities when 
they intersect one of the hexagon's edges and the tendency of all annulus orbits is to have 
p r > at the bisectors. However, it occurs because the Hamiltonian given by (35) and 
(38) is not invariant under the discrete symmetry p t — * — pi, but is only invariant under the 
weaker discrete symmetry (pi,e) — > {~Pi, — e). The parameter e = ±1 is a discrete constant 
of integration that is a measure of the flow of time of the gravitational field relative to the 
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Figure 9: The Poincare plot of the system at H = 1.2 and A = 0. The upper inset shows 
a close up of the upper chaotic region where the Poincare plot is filled in. The lower inset 
shows a close up of the structure in the pretzel region. 
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Figure 10: The Poincare plot for H = 1.2 and A = —0.175, just above the critical value of 
-0.18 for that energy. The upper inset gives a close up of the very tiny remaining region of 
chaos. The lower inset gives a close up of the structure of the lower right pretzel region. 
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Figure 11: The Poincare plot for H = 1.2 and A = 0.06, near the highest value we can 
obtain numerically. The upper and lower insets provide successive closeups of the lower 
right pretzel region showing both the regions of pretzel chaos and the self-similar structure 
at increasingly small scales. 
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particle momenta. In our study we have chosen e = +1 throughout, which has the effect of 
making the principal features of the Poincare plot "squashed" towards the lower right-hand 
side of the figure. If we had chosen e = — 1 this deformation would be towards the lower 
left side. This is reminiscent of the situation for two particles, in which the gravitational 
coupling to the kinetic energy of the particles causes a distorition of the trajectory from an 
otherwise symmetric pattern [19, 20] becoming more pronounced as H increases. 

The general structure of the Poincare section can be seen in figure (9). There is a fixed 
point just right of the center of the plot surrounded by triangular rings. These near-integrable 
curves correspond to the annulus orbits for that energy. All of the possible annulus orbits 
are contained in the largest triangular ring surrounding this region. Its boundary contains a 
thin region of chaos and outside of that, at the bottom and top corners, we have the regions 
corresponding to the pretzel orbits. 

The structures in the lower part and upper corners of the Poincare section is extremely 
complicated. We find a self similarity in these regions with a series of circles bounding the 
quasiperiodic near-integrable regions repeating themselves on increasingly smaller scales. 
We find that each of these circular patterns corresponds to a family of pretzel orbits. The 
two largest circles just below the annulus region correspond to the boomerang-shaped orbits 



The next set of three circles corresponds to the sequence yA 2 B 3 J , and so on. We see 

a collection of crescents between these sets of circles and they correspond to the sequences 
AB 3 A 2 B 3 , AB 3 AB 3 A 2 B 3 , and so on. Inside these circles there is actually a continuum 
of possible circles with a diameter that depends upon the initial conditions. These circles 
are centered around a single point that corresponds to the periodic orbit that these orbits 
approach. 

As H increases, we remarkably do not find a breakdown from regular to chaotic motion 
in this system. As described before we see the plot skewed farther to the bottom right with 
larger H but the topology of the system seems to remain the same. For instance in the 
lower region of the Poincare map we see the same patterns of series of circles at all different 
energies with no sizeable connected areas of chaos. We see the characteristic regions of chaos 
appearing in the regions between the annulus region and the pretzel regions, but we see little 
change in the relative sizes of the chaotic regions with H . 

When we introduce a negative cosmological constant we see some significant changes in 
the characteristics of the Poincare sections shown in figures (10) and (13). First, we see that 
the asymmetric skewing of the graph to the bottom right is reduced and "pushed" back to the 
center. Specifically as A approaches its negative critical value, the Poincare section becomes 
symmetric about the p r = line. In the bottom pretzel region we see that the repeating 
pattern of circles is continued but again centered on the graph instead of skewed to the 
right. This is similar to the low energy Newtonian case studied in [8] where the Poincare 
section is also centered around p r = 0. This is commensurate with the interpretation that 
the negative critical value of A for a given energy corresponds to the point where all of the 
energy in the system is vacuum energy, leaving no (ie very little) energy left for the motion 
of the particles. 

The most suprising feature of the Poincare section when A < is the rapid disappearence 
of all the regions of chaos. In the region between the annulus and pretzel regions we find that 
once-chaotic orbits split into purely annulus or purely pretzel class orbits as A approaches 
its negative critical value. We find no evidence of the onset of a KAM breakdown anywhere 
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Figure 12: The Poincare plot of the system with H = 1.5 and A = 0. The upper inset 
provides a closeup of the upper chaotic region while the lower inset provides a closeup of the 
lower right pretzel region. 
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Figure 13: The Poincare plot for the system at H = 1.5 and A = —0.25, just above the 
critical value at that energy. The upper inset provides a closeup of the upper chaotic region 
which has essentially dissappeared. The lower inset provides a closeup of the lower right 
pretzel region. 
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Figure 14: The Poincare plot of the system at H = 1.5 and A = 0.02, close to the upper 
limit of A that is numerically obtainable. 
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in the plot. This feature is strikingly different from what we have seen in the corresponding 
Newtonian system where chaos is present at all energies. 

In our study of the A < case we were able to numerically solve the equations of motion 
of the hex-particle for all values of A up to the physical limit of the critical negative value 
where all of the energy of the system is "used up" in the contraction of the spacetime. 
However, for the A > case, due to the reduction of V rc as A increases (as previously 
described) we find that the system that was in a perturbative regime quickly changes into an 
intrinsically nonperturbative system and our numerical simulations break down. This places 
much greater restrictions on what values of A > we can test, especially at higher energies. 
Because of this numerical limit, our study of the A > case is somewhat more limited since 
we cannot follow the changes to this system up to a physical limit. However, we have seen 
the beginnings of a number of characteristic changes to the Poincare section as A increases, 
as shown in figures (11) and (14). 

First, many features are skewed farther to the right as A increases. Furthermore many 
of the features in the bottom pretzel region seem to be skewed up and to the right instead 
of down. Below the two large circles in the pretzel region we again see the repeating pattern 
of circles very well defined. 

We also find that the chaotic region in between the boundary of the annulus and pretzel 
regions is also greatly increased, both in the corners of the boundary and along the sides. 
These characteristics are found at all numerically obtainable values of H, though for larger H, 
as described before we cannot increase A as much before the system becomes non-integrable. 

Furthermore, as A increases we at first see the onset of KAM breakdown in the pretzel 
region, followed by the appearence of major regions of chaos. This transition is illustrated in 
figure (15). For H = 1.2, at A = 0.025 we already begin to see some of the lines widen and 
small regions of chaos appear between the groups of ellipses. In these narrow regions the hex- 
particle seems to switch between orbits with an infinitely repeating symbol sequence (which 
correspond to the groups of ellipses) and orbits which undergo additional A type motions in 
quasiperiodic intervals (corresponding to the wavy lines) at seemingly irregular intervals. It 
is these regions that result in the chaotic pretzel orbits mentioned before. However, at this 
point most of the orbits still show very regular motions. As A increases to 0.05 we see these 
new regions of chaos expand around the groups of ellipses. Once we increase A to 0.06 most 
of the lower region of the Poincare section becomes chaotic, with only a few non-connected 
regions of regular motion remaining. This is in sharp contrast with the behaviour of the 
system as A becomes negative, in which no evidence of any KAM breakdown is apparent. 

Unfortunately, due to our previously mentioned limitations on studying systems with 
large cosmological constants, for systems with higher energies we have not been able to 
numerically investigate systems with as large a value of A, so we have not been able to 
observe this full KAM breakdown at all energies. However, for lower values of A we have 
observed the same initial trends of the lines between the series of ellipses begin to thicken 
as shown in the lower right insert of figure (14). Figure (16) shows another example of these 
trends at higher energies. As A is increased, we again see an increase in the area of the upper 
chaotic regions and the appearance of very small regions of chaos in between the series of 
ellipses. This suggests that the preliminary stages of KAM breakdown occur at these higher 
energies as well. 
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Figure 15: A closeup of the pretzel region of the Poincare plots with H = 1.2 for increasing 
values of A. The diagrams are all of the same part of the section. We see clear evidence of 
KAM breakdown as A increases. 
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Figure 16: Regions of the Poincare section at H = 1.8 and A = for the top graphs and 
A = 0.0085 for the bottom graphs. The graphs on the left show closeups of the upper chaotic 
region and the graphs on the right show closeups of the lower pretzel region. 
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8 Discussion 



Our investigation of the cosmological 3-body problem has revealed a number of interesting 
features. We find in general that the presence of a cosmological constant significantly modifies 
the chaotic properties of the relativistic 3-body system, and this in markedly different ways 
depending on its sign. 

For a negative cosmological constant we find that there is a rapid decrease in the amount 
of chaos for all values of H that we were able to investigate. The size of the chaotic regions in 
the Poincare plot are even smaller than in its non-relativistic counterpart, despite the high 
degree of non-linearity in the cosmological system. Indeed, as the cosmological constant 
approaches its negative critical value (defined in 47), the chaotic regions nearly vanish. We 
conjecture that this occurs for arbitrarily large H, motivated primarily by our observation 
that the area of the chaotic regions in the Poincare section seems to be roughly proportional 



to 



for all energies we were able to numerically investigate. 



A-negcrit 

Conversely, we find an increase in the area of chaotic regions in the Poincare section 
when the cosmological constant is positive, both in the regions between the annulus and 
pretzel orbits, and within the regions corresponding to the pretzel orbits. Though we were 
unable to investigate very large cosmological constants for systems with higher energies, at 
higher energies, even for the reasonably small positive cosmological we were able to study, we 
observed the lines in the pretzel regions of the Poincare section starting to thicken and very 
small regions of chaos appear between the groups of ellipses , reminicent of the preliminary 
stages of KAM breakdown. We therefore conjecture that this increase in chaos occurs for 
positive A at all energies. 

We close with some comments on future work. The unequal mass case remains to be 
explored. While we expect that the general feature of chaos increasing/decreasing with 
positive /negative A will still be present, there could be a number of surprising features in 
the details relative to the A = case [9]. However our primary limitation has been that of 
exploring energies that are below the critical value of the potential. This same limitation 
was present in previous studies of the A = 3-body system [8, 9] - to move beyond it will 
require employing a time parameter that is not coordinate time, as well as more sophisticated 
numerical algorithms that avoid instabilities we encountered at higher energies. 
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